NASA/ CR-20 1 2-2 17562 
NIA Report No. 2012-01 



[National 

I INSTITUTE O F 

Aerospace 



Development and Application of Benchmark 
Examples for Mixed-Mode I/II Quasi- Static 
Delamination Propagation Predictions 

Ronald Krueger 

National Institute of Aerospace, Hampton, Virginia 


April 2012 


NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to 
the advancement of aeronautics and space science. 

The NASA scientific and technical information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI program operates under the 
auspices of the Agency Chief Information Officer. It 
collects, organizes, provides for archiving, and 
disseminates NASA’s STI. The NASA STI program 
provides access to the NASA Aeronautics and Space 
Database and its public interface, the NASA Technical 
Report Server, thus providing one of the largest 
collections of aeronautical and space science STI in 
the world. Results are published in both non-NASA 
channels and by NASA in the NASA STI Report 
Series, which includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or 
theoretical analysis. Includes compilations of 
significant scientific and technical data and 
information deemed to be of continuing 
reference value. NASA counterpart of peer- 
reviewed formal professional papers, but having 
less stringent limitations on manuscript length 
and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or of 
specialized interest, e.g., quick release reports, 
working papers, and bibliographies that contain 
minimal annotation. Does not contain extensive 
analysis. 

• CONTRACTOR REPORT. Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from NASA 
programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’s mission. 

Specialized services also include creating custom 

thesauri, building customized databases, and 

organizing and publishing research results. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page at 
http://www.sti. nasa.gov 

• E-mail your question via the Internet to 
help@sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 
at 443-757-5803 

• Phone the NASA STI Help Desk at 
443_757_5802 

• Write to: 

NASA STI Help Desk 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 


NASA/ CR-20 1 2-2 17562 
NIA Report No. 2012-01 



National 
Institute of 
Aerospace 



Development and Application of Benchmark 
Examples for Mixed-Mode I/II Quasi- Static 
Delamination Propagation Predictions 

Ronald Krueger 

National Institute of Aerospace, Hampton, Virginia 


National Aeronautics and 
Space Administration 

Prepared for Langley Research Center 
Under Cooperative Agreement NNL09AA00A 


Langley Research Center 
Hampton, Virginia 23681-2199 


April 2012 


The use of trademarks or names of manufacturers in this report is for accurate reporting and does not constitute an 
official endorsement, either expressed or implied, of such products or manufacturers by the National Aeronautics 
and Space Administration. 


Available from: 

NASA Center for AeroSpace Information 
7115 Standard Drive 
Hanover, MD 21076-1320 
443-757-5802 


DEVELOPMENT AND APPLICATION OF BENCHMARK EXAMPLES FOR MIXED- 
MODE I/II QUASI-STATIC DELAMINATION PROPAGATION PREDICTIONS 


Ronald Krueger* 

ABSTRACT 

The development of benchmark examples for quasi-static delamination propagation 
prediction is presented and demonstrated for a commercial code. The examples are based on 
finite element models of the Mixed-Mode Bending (MMB) specimen. The examples are 
independent of the analysis software used and allow the assessment of the automated 
delamination propagation prediction capability in commercial finite element codes based on 
the virtual crack closure technique (VCCT). First, quasi-static benchmark examples were 
created for the specimen. Second, starting from an initially straight front, the delamination 
was allowed to propagate under quasi-static loading. Third, the load-displacement 
relationship from a propagation analysis and the benchmark results were compared, and 
good agreement could be achieved by selecting the appropriate input parameters. Good 
agreement between the results obtained from the automated propagation analysis and the 
benchmark results could be achieved by selecting input parameters that had previously 
been determined during analyses of mode I Double Cantilever Beam and mode II End 
Notched Flexure specimens. The benchmarking procedure proved valuable by 
highlighting the issues associated with choosing the input parameters of the particular 
implementation. Overall, the results are encouraging, but further assessment for mixed- 
mode delamination fatigue onset and growth is required. 


1. INTRODUCTION 

Over the past two decades, the use of fracture mechanics has become common practice to 
characterize the onset and growth of delaminations. In order to predict delamination onset or 
growth, the calculated strain energy release rate components are compared to interlaminar fracture 
toughness properties measured over a range from pure mode I loading to pure mode II loading. 

The virtual crack closure technique (VCCT) is widely used for computing energy release rates 
based on results from continuum (2D) and solid (3D) finite element (FE) analyses and to supply the 
mode separation required when using the mixed-mode fracture criterion [1, 2], The virtual crack 
closure technique was recently implemented into several commercial finite element codes. As new 
methods for analyzing composite delamination are incorporated into finite element codes, the need 
for comparison and benchmarking becomes important since each code requires specific input 
parameters unique to its implementation. These parameters are unique to the numerical approach 
chosen and do not reflect real physical differences in delamination behavior. 

An approach for assessing the mode I, mode II and mixed-mode I and II, delamination 
propagation capabilities in commercial finite element codes under static loading was recently 
presented and demonstrated for VCCT for ABAQUS® * 1 [3,4] as well as MD Nastran™ and Marc™ 2 
[5]. First, benchmark results were created manually for finite element models of the mode I Double 
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Cantilever Beam (DCB), the mode II End Notched Flexure (ENF) and the mixed-mode I/II Single 
Leg Bending (SLB) specimen. Second, starting from an initially straight front, the delamination was 
allowed to propagate using the automated procedure implemented in the finite element software. 
The approach was then extended to allow the assessment of the delamination fatigue growth 
prediction capabilities in commercial finite element codes [4,6]. As for the static case, benchmark 
results were created manually first for the mode I Double Cantilever Beam (DCB) and the mode II 
End Notched Flexure (ENF) specimen. Second, the delamination was allowed to grow under cyclic 
loading in a finite element model of a commercial code. In general, good agreement between the 
results obtained from the propagation and growth analysis and the benchmark results could be 
achieved by selecting the appropriate input parameters. Overall, the results were encouraging but 
showed that additional assessment for mixed-mode delamination cases is required. 

The objective of the present study was to create additional benchmark examples, independent of 
the analysis software used, which allows the assessment of the quasi-static delamination 
propagation prediction capabilities in commercial finite element codes. For the simulation of mixed- 
mode 1/11 fracture, the Mixed-Mode Bending (MMB) specimen was selected as shown in Figure 1. 
Dimensions, layup and material properties were taken directly from a related experimental study 

Ul 

Static benchmark results were created for three mixed-mode ratios {Gii/Gt =0.2, 0,5 and 0.8) 
based on the approach developed earlier [3]. To create the benchmark results, two-dimensional 
finite element models were used for simulating the MMB specimens with different delamination 
lengths ao. For each delamination length modeled, the load, P, and the displacement, w, were 
monitored. The total strain energy release rate, Gr, and mixed-mode ratio, Gii/Gt, were calculated 
for a fixed applied displacement. It is assumed that the delamination propagates when the computed 
energy release rate, Gt, reaches the mixed-mode fracture toughness G c . Thus, critical loads and 
critical displacements for delamination propagation were calculated for each delamination length 
modeled. From these critical load/displacement results, benchmark solutions were created. It is 
assumed that the load/displacement relationship computed during automatic propagation should 
closely match the benchmark cases. 

After creating the benchmark cases, the approach was demonstrated for the commercial finite 
element code ABAQUS®. Starting from an initially straight front, the delamination was allowed to 
propagate under quasi-static loading based on the algorithms implemented into the software. Input 
control parameters were varied to study the effect on the computed delamination propagation. The 
benchmark enabled the selection of the appropriate input parameters that yielded good agreement 
between the results obtained from the propagation analysis and the benchmark results. Once the 
parameters have been identified, they may then be used with confidence to model delamination 
growth for more complex configurations. 

In this paper, the development of the benchmark cases for the assessment of the quasi-static 
delamination propagation is presented. Examples of automated propagation analyses are shown, and 
the selection of the required code specific input parameters are discussed. 


2. METHODOLOGY BASED ON FRACTURE MECHANICS 

For the current numerical investigation, the Mixed-Mode Bending (MMB) specimen, as shown 
in Figure 1, was chosen since it is simple, and exhibits the mixed-mode 1/11 opening fracture mode 
over a wide range of mixed-mode ratios Gii/Gt. For the current investigation, finite element models 
were created for three distinct mixed-mode ratios ( Gii/Gt =0.2, 0,5 and 0.8). The methodology for 
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delamination propagation, was applied to the MMB specimen to create the benchmark example [ 8 , 
9]. For the current study, MMB specimens made of IM7/8552 graphite/epoxy with a unidirectional 
layup, [0]24, were modeled. The material, layup, overall specimen dimensions including initial crack 
length, ao, were identical to specimens used in related experimental studies [7,10,1 1]. The material 
properties are given in Tables I and II. 

A quasi-static mixed-mode fracture criterion is determined by plotting the interlaminar fracture 
toughness, G c , versus the mixed-mode ratio, Gji/Gt as shown in Figure 2 . The fracture criterion is 
generated experimentally using pure Mode I (Gn/Gj=0) Double Cantilever Beam (DCB) tests [11], 
pure Mode II (G//G 7 =l) End-Notched Flexure (ENF) tests [10], and Mixed Mode Bending (MMB) 
tests of varying ratios of G/ and Gn [7]. For the material used in this study, the experimental data 
(open symbols) and mean values (filled symbols) are shown in Figure 2. A 2D fracture criterion was 
suggested by Benzeggah and Kenane [ 12 ] using a simple mathematical relationship between G c and 
Gjj/Gt 


a r -G, c+ (G„ r -G„) 


Ig 

11 

\G T ] 
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In this expression, G/ c and G// c are the experimentally determined fracture toughness data for mode I 
and II as shown in Figure 2. The factor 17 was determined by a curve fit using the Levenberg- 
Marquardt algorithm in the KaleidaGraph™ graphing and data analysis software [13]. The 
parameters G/ c , G// c and 17 are required input to perform a VCCT analysis in ABAQUS® as 
discussed in the appendix. 


3. FINITE ELEMENT MODELING 
3.1 Model description 

An example of a two-dimensional finite element model of a Mixed-Mode Bending (MMB) 
specimen with boundary conditions is shown in Figure 3 for a mixed-mode ratio G///G 7 - =0.2. Based 
on previous experience [3,4,6], the specimen was modeled with solid plane strain elements (CPE4I) 
in ABAQUS® Standard 6.10. Along the length, all models were divided into different sections with 
different mesh refinement as shown in Figure 3 a. The MMB specimen was modeled with six 
elements through the specimen thickness (2h) as shown in the detail of Figure 3b. The resulting 
element length at the delamination tip was Aa=0.5 mm. The load apparatus was modeled explicitly 
using rigid beam elements (R2D2) as shown in Figure 3 a. Multi-point constraints were used to 
connect the rigid elements with the planar model of the specimen and enforce the appropriate 
boundary conditions as shown in Figures 3b and c. 

The plane of delamination was modeled as a discrete discontinuity in the center of the specimen. 
For the analysis with ABAQUS® 6.10, the models were created as separate meshes for the upper 
and lower part of the specimens with identical nodal point coordinates in the plane of delamination 
[14], Two surfaces (top and bottom surface) were defined to identify the contact area in the plane of 
delamination as shown in Figures 3b and c. Additionally, a node set was created to define the intact 
(bonded nodes) region. Two models, used to simulate mixed-mode ratios Gji/Gt =0.5 and 0 . 8 , are 
shown in Figure 4. The mesh of the specimen was kept the same for all three mode ratios. Only the 
lengths of the rigid elements used to simulate the load apparatus were changed as shown in 
Figure 4. 
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Examples of three-dimensional finite element models of the MMB specimen are shown in 
Figures 5 and 6. Along the length, all models were divided into different sections with different 
mesh refinement. A refined mesh was used in the center of the MMB specimen as shown in the 
detail of Figure 5b. Across the width, a uniform mesh (25 elements) was used to avoid potential 
problems at the transition between a coarse and finer mesh [3-6]. Through the specimen thickness 
(2/z), six elements were used as shown in the detail of Figure 5b. The resulting element length at the 
delamination tip was Aa=0.5 mm. The specimen was modeled with solid brick elements (C3D8I), 
which had yielded excellent results in previous studies [3,4,6]. The load apparatus was modeled 
explicitly using rigid plate elements (R3D4) as shown in Figure 5a. As before for the two- 
dimensional model, multi-point constraints were used to connect the rigid elements with the solid 
model of the specimen and enforce the appropriate boundary conditions. Two models used to 
simulate mixed-mode ratios Gi/Gt =0.5 and 0.8 are shown in Figure 6. The mesh of the specimen 
was kept the same for all three mode ratios. Only the lengths of the rigid plate elements used to 
simulate the load apparatus were changed as shown in Figure 6. 

3.2 Static delamination propagation analysis 

For the automated delamination propagation analysis, the YCCT implementation in ABAQFFS® 
Standard 6.10 was used. The plane of delamination in three-dimensional analyses is modeled using 
the existing AB AQU S ®/S tan dard crack propagation capability based on the contact pair capability 
[14]. Additional element definitions are not required, and the underlying finite element mesh and 
model does not have to be modified [14]. The implementation offers a crack and delamination 
propagation capability in ABAQFFS®. It is implied that the energy release rate at the crack tip is 
calculated at the end of a converged increment. Once the energy release rate exceeds the critical 
strain energy release rate (including the user-specified mixed-mode criteria as shown in Figure 2), 
the node at the crack tip is released in the following increment, which allows the crack to propagate. 
To avoid sudden loss of stability when the crack tip is propagated, the force at the crack tip before 
advancement is released gradually during succeeding increments in such a way that the force is 
brought to zero no later than the time at which the next node along the crack path begins to open 
[14,15]. 

In addition to the mixed-mode fracture criterion, YCCT for ABAQFFS® requires additional input 
for the propagation analysis. If a user specified release tolerance is exceeded in an increment 
(G-G c )/G c > release tolerance, a cutback operation is performed which reduces the time increment. 
In the new smaller increment, the strain energy release rates are recalculated and compared to the 
user specified release tolerance. The cutback reduces the degree of overshoot and improves the 
accuracy of the local solution [14]. A release tolerance of 0.2 is suggested in the handbook [14]. To 
help overcome convergence issues during the propagation analysis, ABAQFFS® provides: 

• contact stabilization which is applied across only selected contact pairs and used to 
control the motion of two contact pairs while they approach each other in multi-body 
contact. The damping is applied when bonded contact pairs debond and move away 
from each other [14,15] 

• automatic or static stabilization which is applied to the motion of the entire model and is 
commonly used in models that exhibit statically unstable behavior such as buckling 
[14,15] 

• viscous regularization, which is applied only to nodes on contact pairs that have just 
debonded. The viscous regularization damping causes the tangent stiffness matrix of the 
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softening material to be positive for sufficiently small time increments. Viscous 
regularization damping in VCCT for ABAQUS® is similar to the viscous regularization 
damping provided for cohesive elements and the concrete material model in 
ABAQUS®/Standard [14,15]. Further details about the required input parameters are 
discussed in the appendix. 

For automated propagation analysis, it was assumed that the computed behavior should closely 
match the benchmark results created below. For all analyses, the elastic constants (given in Table I) 
and the input to define the fracture criterion (given in Table II) were kept constant. Based on 
previous results [3], automatic or static stabilization was not used in this study. To limit the 
scope of the study, viscous regularization was also not considered. Instead, the study focused on 
using the default values and parameter combinations that had recently yielded good results [4], 
This approach was used to gain confidence in the parameters that had been identified previously 
and thus determine if the selection of the parameters could be considered problem independent. 
Therefore, only the following items were varied to study the effect on the automated delamination 
propagation behavior during the analysis. 

• The release tolerance ( reltol) was varied. 

• Analyses were performed with and without contact stabilization. For analyses that 
included contact stabilization only, a single stabilization factor (c.s= 1 x 1 0“ 6 ) was used. 

• Two- and three-dimensional models with different types of elements were used. 


4. DEVELOPMENT OF THE STATIC BENCHMARK CASES 
4.1 Benchmark case for 20% mode II 

The static benchmark case was created based on the approach developed earlier [3]. Two- 
dimensional finite element models simulating MMB specimens with 17 different delamination 
lengths ao were created (25.4 mm< ao< 70.6 mm). For each delamination length modeled, the load, 
P, and displacement, vv, were monitored as shown in Figure 7 (colored lines) for the case of 20% 
mode II ( Gi/Gt =0.2). Using VCCT, the total strain energy release rate, Gt, and the mixed-mode 
ratio Gji/Gt were computed at the end of the analysis as shown in Figure 7. The failure index Gt/G c 
was calculated by correlating the computed total energy release rate, Gt, with the mixed-mode 
fracture toughness, G c , of the graphite/epoxy material. As discussed before (see Figure 2), the 
mixed-mode fracture toughness, G c , is a function of the mixed-mode ratio Gji/Gt. Hence, the mixed- 
mode fracture toughness, G c for each computed mixed-mode ratio ( Gji/Gt —0.2) was obtained from 
the curve fit of the material data (solid red curve) as illustrated in Figure 8 (blue arrows). It is 
assumed that the delamination propagates when the failure index Gt/G c reaches unity. Therefore, 
the critical load, P crit , can be calculated based on the relationship between load, P, and the energy 
release rate, G [16], 


p 2 dc p 
2 dA 


( 2 ) 


In equation (2), Cp is the compliance of the specimen, and d A is the increase in surface area 
corresponding to an incremental increase in load or displacement at fracture. The critical load, P cri t, 
and critical displacement, w crit , were calculated for each delamination length modeled 
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and the results were included in the load/displacement plots as shown in Figure 9 (solid red circles). 
These critical load/displacement results indicated that, with increasing delamination length, less 
load is required to extend the delamination up to ao ~48 mm. To further extend the delamination an 
increase in load is required past w= 3.5 mm. 

As shown in Figure 10, these critical load/displacement results (solid black dots) can be used to 
create a benchmark solution (solid red line) for applied displacements, w. The benchmark result 
may also be visualized by plotting the prescribed displacements, w, at delamination growth onset 
versus the increase in delamination length, a*, as illustrated in Figure 11. It is assumed that the 
load/displacement, or applied displacement/delamination-length relationship computed during 
automatic propagation should closely match the benchmark cases. 

4.2 Benchmark case for 50% mode II 

The benchmark case for 50% mode II ( Gi/Gt =0.5) was created similarly. Two-dimensional 
finite element models simulating MMB specimens with 16 different delamination lengths ao were 
created (25.4 mm< ao< 70.6 mm) to study the case of Gi/Gt =0.5. For each delamination length 
modeled, the load, P, and displacement, w, were monitored as shown in Figure 12 (colored lines). 
Using YCCT, the total energy release rate, Gt, and the mixed-mode ratio Gi/Gt were computed at 
the end of the analysis as shown in Figure 12. The mixed-mode fracture toughness, G c for each 
computed mixed-mode ratio ( Gi/Gt ~0.5) was obtained from the curve fit of the material data (solid 
red curve) as illustrated previously in Figure 8 (green arrows). The critical load, P crit , and critical 
displacement, w crih were calculated for each delamination length modeled using equation (3), and 
the results were included in the load/displacement plots as shown in Figure 13 (solid red circles). 
These critical load/displacement results indicated that, with increasing delamination length, less 
load is required to extend the delamination up to ao ~48 mm. To further extend the delamination an 
increase in load is required past w= 2.0 mm. For the first three delamination lengths, ao, plotted in 
Figure 13, the values of the critical displacements also decreased at the same time. This means that 
the MMB specimen exhibits unstable delamination propagation under load control as well as 
displacement control in this region for Gi/Gt =0.5. The remaining critical load/displacement results 
indicated stable propagation. From these critical load/displacement results (solid grey circles and 
dashed line), two benchmark solutions can be created as shown in Figure 14. During the analysis, 
either prescribed displacements, w, or nodal point loads, P, are applied. For the case of prescribed 
displacements, w, (dashed blue line), the applied displacement must be held constant over several 
increments once the critical point (Peru, Wa-u) is reached, and the delamination front is advanced 
during these increments. Once the critical path (dashed grey line) is reached, the applied 
displacement is increased again incrementally. For the case of applied nodal point loads (dashed red 
line), the applied load must be held constant while the delamination front is advanced during these 
increments. Once the critical path (dashed grey line) is reached, the applied load is increased again 
incrementally. 

The benchmark result for prescribed displacements may also be visualized by plotting the 
prescribed displacements, w, at delamination growth onset versus the increase in delamination 
length, a*, as illustrated in Figure 15. For the case of applied nodal point loads, the benchmark 
result may be visualized by plotting the applied loads, P, at delamination growth onset versus the 
increase in delamination length, a* as illustrated in Figure 16. It is assumed that the 
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load/displacement, load/delamination-length or applied displacement/delamination-length 
relationship computed during automatic propagation should closely match the benchmark cases. 

4.3 Benchmark case for 80% mode II 

The benchmark case for 80% mode II ( Gu/Gt =0.8) was created based on the approach 
developed earlier [3] and discussed above for 20% and 50% mode II. Two-dimensional finite 
element models simulating MMB specimens with 15 different delamination lengths ao were created 
(25.4 mm< ao<10.6 mm) to study the case of Gi/Gt =0.8. For each delamination length modeled, 
the load, P, and displacement, w, were monitored as shown in Figure 17 (colored lines). Using 
VCCT, the total energy release rate, GV, and the mixed-mode ratio Gu/Gt were computed at the end 
of the analysis as shown in Figure 17. The mixed-mode fracture toughness, G c for each computed 
mixed-mode ratio ( Gu/Gt —0.8) was obtained from the curve fit of the material data (solid red 
curve) as illustrated previously in Figure 8 (purple arrows). The critical load, P crit , and critical 
displacement, w crih were calculated for each delamination length modeled, using equation (3), and 
the results were included in the load/displacement plots as shown in Figure 18 (solid red circles). 

These critical load/displacement results indicated that, with increasing delamination length, less 
load is required to extend the delamination up to ao ~48 mm. To further extend the delamination an 
increase in load is required past vv=2.0 mm.For the first five delamination lengths, ao, plotted in 
Figure 18, the values of the critical displacements also decreased at the same time. This means that 
the MMB specimen exhibits unstable delamination propagation under load control as well as 
displacement control in this region for Gu/Gt =0.8. The remaining critical load/displacement results 
indicated stable propagation. From these critical load/displacement results (solid grey circles and 
dashed line), two benchmark solutions can be created as shown in Figure 19. During the analysis, 
either prescribed displacements, w, or nodal point loads, P, are applied. For the case of prescribed 
displacements, w, (dashed blue line), the applied displacement must be held constant over several 
increments once the critical point (Pent, Wa-n) is reached, and the delamination front is advanced 
during these increments. Once the critical path (dashed grey line) is reached, the applied 
displacement is increased again incrementally. For the case of applied nodal point loads (dashed red 
line), the applied load must be held constant while the delamination front is advanced during these 
increments. Once the critical path (dashed grey line) is reached, the applied load is increased again 
incrementally. 

The benchmark result for prescribed displacements may also be visualized by plotting the 
prescribed displacements, w, at delamination growth onset versus the increase in delamination 
length, a*, as illustrated in Figure 20. For the case of applied nodal point loads, the benchmark 
result may be visualized by plotting the applied loads, P, at delamination growth onset versus the 
increase in delamination length, a* as illustrated in Figure 21. It is assumed that the 
load/displacement, load/delamination-length or applied displacement/delamination-length 
relationship computed during automatic propagation should closely match the benchmark cases. 


5. STATIC ANALYSIS BENCHMARKING 

5.1 Computed energy release rates and mixed-mode ratios 

First, static analyses were performed with propagation disabled to compute the total energy 
release rate, Gt, the mode II component, Gu, and the mixed-mode ratio Gu/Gt. For this static 
analysis, an applied maximum displacement of w=2.0 mm was chosen for all models. To ensure that 
the model geometry and material input data for all models produced consistent results, the 
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computed results were evaluated. The results were obtained from models shown in Figures 3 
through 6. 

5.1.1 Comparison of results obtained from two- and three-dimensional models 

For all mixed-mode ratios (results for Gu/Gi=()2 in blue; 0.5 in green and 0.8 in black) the two- 
dimensional models (results with open symbols and solid lines) yield total energy release rates, Gt, 
that are almost identical to the results obtained from three-dimensional models (solid symbols and 
dashed lines) as shown in Figure 22. Comparing the computed mode II strain energy release rates, 
Gjj (plotted in Figure 23), and the mixed-mode ratios, Gu/Gt, (plotted in Figure 24), confirms the 
good agreement between the different models. The computed mixed-mode ratios, Gu/Gt, are 
slightly lower than the targeted values for the test (dashed lines) [7] a trend that appears to increase 
with increasing applied displacement, w. Based on the current findings, it was assumed that results 
obtained from two- and three-dimensional models computed later during the assessment phase were 
easily comparable. 

For each of the two-dimensional finite element models representing different delamination 
lengths (25.4 mm< ao<162 mm), the mixed-mode ratios corresponding to an applied displacement 
w=2.0 mm, were computed. The results were plotted versus the increase in delamination length, a* 
as illustrated in Figure 25. For all three mixed-mode ratios studied (results for Gu/Gt=02 in blue; 
0.5 in green and 0.8 in black) the computed mixed-mode ratios increase slightly with increasing 
delamination length until a*~ 20 mm when the mixed-mode ratios reach the targeted values for the 
test (dashed lines) [7]. The computed mixed-mode ratios increased further until the delamination tip 
is located under the loading point. For delaminations extending beyond the load point, the computed 
mixed-mode ratios increased for models with target values Gu/Gr=0.5 and 0.8 and decreased for 
models with target values Gu/Gr=02. It should be noted, that delaminations extending beyond the 
load point could be included in the current benchmarking process, which is entirely geared toward 
numerical analysis. For actual physical tests geared towards obtaining mixed-mode fracture 
toughness values however, delaminations extending beyond the load point become meaningless. 

5.1.2 Computed distributions across the specimen width 

For all the mixed-mode ratios studied (Gu/Gf=02, 0.5 and 0.8), the total energy release rate, Gt, 
the mode II component, Gu, the mixed-mode ratio, Gu/Gt, and the failure index, Gt/G c , were also 
computed and plotted versus the normalized width, y/B, as shown in Figures 26 through 29. For a 
delamination length ao= 25.4 mm the computed total energy release rates rates, Gt (results for 
G i //G j=0 2 blue circles; 0.5 green squares and 0.8 black diamonds), were plotted versus the 
normalized width, y/B, of the specimen as shown in Figure 26. The results obtained from two- 
dimensional finite element models (open symbols) are plotted at the center of the specimen y/B =0. 
The values obtained from two-dimensional and three-dimensional analyses are in good agreement 
as discussed above (see results plotted in Figures 22 through 25). The results obtained from three- 
dimensional models (solid symbols) indicate that qualitatively, the total strain energy release rate is 
fairly constant over almost the entire width of the specimen. Only near the edges of the specimen is 
a decrease in total energy release rate observed for the cases Gu/Gt=02 and 0.5 and an increase in 
total energy release rate is observed for the mode II dominated case Gu/Gt=0.S. The mode II 
components, Gu, for each case were plotted in Figure 27. An increase in Gu is visible near the edges 
for the mode II dominated case Gu/Gj=0.S and also for the equal mode case Gu/Gt=0.5. The mixed- 
mode ratios, Gu/Gt, plotted in Figure 28 indicate a relative dominance of Gu near the edges of the 
specimens compared to the distribution in the center of the specimens. For all three-dimensional 
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models, the failure indices, Gt/G c , plotted in Figure 29, were somewhat higher across the entire 
width than the results from two-dimensional planar finite element models. The failure indices 
dropped near the edges, which suggested that delamination propagation was expected to locally lag 
behind and turn the initially straight front into a thumbnail front for Gji/Gt =0.2 and 0.5 similar to 
fronts observed in a Double Cantilever Beam (DCB) specimen. For Gi/Gt =0.8, the drop of the 
failure index near the edges was less pronounced and a straight front across the entire width was 
expected. The computed front shapes will be discussed in detail in a later section. 

5.2 Results from automated delamination propagation analysis for 20% mode II 

The propagation analysis was performed in two steps using the models shown in Figures 3 and 5 
starting from an initial delamination length, ao= 25.4 mm. In the first step, a displacement of w=l .5 
mm was applied which nearly equaled the critical displacement, w crit =\ .64 mm, determined earlier. 
In the second step, the applied displacement was increased to w=8.0 mm. For this second step, 
automatic incrementation was used in ABAQUS® and a small increment size (0.5% of the total 
step) was chosen at the beginning of the step. To reduce the risk of numerical instability and early 
termination of the analysis, a minimum allowed increment size (10‘ of the total step) was also 
chosen. The analysis was limited to 5000 increments. 

Initially, analyses were performed using two-dimensional planar models (shown in Figure 3) 
without stabilization or viscous regularization. Based on previous results [3,4], release tolerance 
values smaller than the default value (reltol= 0.2) [14] were used since these values yielded better 
results. Using reltol= 0.1, the load dropped at the critical point, however the computed path stayed 
above the benchmark as shown in Figure 30 (solid blue line) where the computed resultant force 
(load P ) is plotted versus the applied displacement w. As the displacement continued to increase 
over 3.0 mm, the computed load/displacement path converged to the benchmark result (solid grey 
circles and solid grey line). To improve the results, the release tolerance was decreased to 
reltol= 0.01. The results (solid red line) matched the benchmark case. Due to the fine mesh, only a 
small saw-tooth pattern was observed. 

Based on problems identified during previous analyses, automated or static stabilization was not 
used in this study [3]. The results computed when contact stabilization (cs) was added are plotted in 
Figure 31. For a small stabilization factor (cs=lxl 0" 6 ) and a release tolerance suggested in the 
handbook (reltol= 0.2) [14], the load decreased, and delamination propagation started shortly after 
reaching the critical point of the benchmark solution (solid green line). The load/displacement path 
then ran parallel to the constant displacement branch of the benchmark result and later converged to 
the benchmark result. To reduce the observed overshoot, the release tolerance was reduced. For a 
stabilization factor of cv=lxlO" 6 and release tolerance reltol= 0.1 (solid blue line), the overshoot was 
reduced, and the computed load/displacement path then ran closer to the constant displacement 
branch of the benchmark result. This result was practically identical to the results discussed above 
without stabilization. As above, further reducing the release tolerance (reltol= 0.01) yielded results 
that were in excellent agreement with the benchmark (solid red line) and due to the fine mesh, only 
a small saw-tooth pattern was observed. 

An alternative way to plot the benchmark is shown in Figure 32 where the applied displacement 
w is plotted versus the increase in delamination length a*. This way of presenting the results is 
shown because it may be of advantage for large structures where local delamination propagation 
may have little effect on the global stiffness of the structure and may therefore not be visible in a 
global load/displacement plot. However, extracting the delamination length a from the finite 
element results required more manual, time consuming post-processing of the results compared to 
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the relatively simple and readily available output of nodal displacements and forces. The results 
plotted in Figure 30 are the examples that were discussed above and were shown in the global 
load/displacement plots of Figure 31. The conclusions that can be drawn from this plot are identical 
to those discussed above. 

Results obtained from three-dimensional models are shown in Figures 33 and 34. Based on the 
results from two-dimensional planar models shown previously, contact stabilization was added to 
the analyses to help overcome convergence issues. The results (solid blue line) computed for a 
release tolerance value (reltol=(). 1 ) and a small stabilization factor (cs=lxl0' 6 ) are plotted in 
Figure 33. The load dropped and delamination propagation started upon reaching the critical point 
of the benchmark solution (solid circles and solid grey line). The load/displacement path then ran 
parallel to the benchmark result but was shifted slightly towards lower loads. Deviation from the 
benchmark may be explained by the fact that the benchmark results were created using two- 
dimensional planar finite element models of the MMB specimen. As shown in Figure 29, the three- 
dimensional models yield a failure index which is somewhat higher across the entire width of the 
specimen compared to the results from two-dimensional planar finite element models. Therefore, 
delamination propagation is expected to start prior to the benchmark results obtained from two- 
dimensional planar finite element analysis, which ultimately leads to a shift of the entire results plot 
towards lower loads. 

For three distinct locations across the width of the specimen (two edges and center as shown 
later in Figure 53) the results from above were plotted in Figure 34, where the applied displacement 
w is plotted versus the increase in delamination length a *. The results for all three locations (solid 
green line for edge 1, solid blue line for center location and solid red line for edge 2) are in excellent 
agreement with the benchmark case (solid circles and solid grey line). Any deviation from the 
benchmark results mentioned above is less visible compared to the load/displacement plot of 
Figure 33. 

5.3 Results from automated delamination propagation analysis for 50% mode II 

5.3.1 Computed delamination propagation for applied displacement 

As discussed above, the propagation analysis was performed in two steps using the models 
shown in Figures 4a and 6a starting from an initial delamination length, ao=25.4 mm. In the first 
step, a displacement of w=\2 mm was applied which nearly equaled the critical displacement, 
w cr u= 1.34 mm, determined earlier. In the second step, the applied displacement was increased to 
w=8.0 mm. For this second step, automatic incrementation was used in ABAQUS® and a small 

increment size (0.5% of the total step) was chosen at the beginning of the step. To reduce the risk of 
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numerical instability and early termination of the analysis, a minimum allowed increment size (10' 
of the total step) was also chosen. The analysis was limited to 5000 increments. 

Initially, analyses were performed using two-dimensional planar models (shown in Figure 4a) 
without stabilization or viscous regularization. Release tolerance values equal to the default value 
(reltol={)2) [14] and smaller were used since these values had yielded better results in the past [3,4], 
Using the default value (reltol= 0.2), the computed load and displacement exceeded the critical point 
before the initial load drop occurred and delamination propagation started (solid green line) as 
shown in Figure 35 where the computed resultant force (load P) is plotted versus the applied 
displacement w. The computed path initially stayed above the benchmark, however, as the 
displacement continued to increase over 2.5 mm, the computed load/displacement path converged 
to the benchmark result (solid grey circles and dashed grey line). To reduce the observed overshoot, 
the release tolerance was decreased. For a release tolerance reltol= 0.1, the overshoot was reduced 
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and the results improved (solid blue line) and shifted towards the benchmark case. For a release 
tolerance reltol= 0.01, the analysis terminated once the critical point was reached (solid red line) due 
to convergence problems. 

The results computed when contact stabilization (cs) was added are plotted in Figure 36. Using 
a stabilization factor of c.s=lxl0' 6 , and a release tolerance (reltol= 0.01) yielded results that were in 
excellent agreement with the benchmark (solid red line). Due to the fine mesh, only a small saw- 
tooth pattern was observed. The result confirms previous observations where small release tolerance 
values (reltol= 0.01) in conjunction with a stabilization factor of cs=lxl 0‘ 6 had yielded excellent 
agreement with the benchmark results [3,4], 

An alternative way to plot the benchmark is shown in Figure 37 where the applied displacement 
w is plotted versus the increase in delamination length a*. The results plotted in Figure 37 are the 
examples that were discussed above and were shown in the global load/displacement plots of 
Figures 35 and 36. The conclusions that can be drawn from this plot are identical to those discussed 
above. 

The result obtained from a three-dimensional analysis is shown in Figure 38. Based on the 
results from two-dimensional planar models shown above, contact stabilization was added to the 
analysis to help overcome convergence issues. Since the computer time increases with tighter 
release tolerances, a tolerance value (reltol= 0.1) was selected. The result (solid blue line) computed 
for a release tolerance value reltol= 0.1 and a small stabilization factor cs=lxl0" 6 are plotted in 
Figure 38. The load dropped and delamination propagation started upon reaching the critical point 
of the benchmark solution (solid circles and dashed grey line). The load/displacement path then ran 
parallel to the benchmark result but was shifted slightly towards lower loads. Deviation from the 
benchmark may be explained by the fact that the benchmark results were created using two- 
dimensional planar finite element models of the MMB specimen. As shown in Figure 29, the three- 
dimensional models yield a failure index which is somewhat higher across the entire width than the 
results from two-dimensional planar finite element models. Therefore, delamination propagation is 
expected to start prior to the benchmark results obtained from two-dimensional planar finite element 
analysis, which ultimately leads to a shift of the entire results plot towards lower loads. 

5.3.2 Computed delamination propagation for applied quasi-static load 

The propagation analysis was performed in two steps using the models shown in Figures 4a and 
6a starting from an initial delamination length, ao=25A mm. In the first step, a load P=360 N was 
applied which equaled nearly the critical load, P C nt = 385 N, determined earlier. In the second step, 
the total load was increased (P=600 N). Automatic incrementation was used with a small increment 
size at the beginning (0.5% of the total step) and a very small minimum allowed increment (10' of 
the total step) to reduce the risk of numerical instability and early termination of the analysis. The 
analysis was limited to 5000 increments. 

The same steps discussed in the section on applied displacement were followed. Initially, 
analyses were performed using two-dimensional planar models without stabilization or viscous 
regularization. For the default value reltol= 0.2, the load increase stopped after reaching the critical 
point, but the analysis terminated immediately (solid green line) due to convergence problems as 
shown in Figure 39. The release tolerance was not increased as suggested by the ABAQUS® error in 
the message (.msg) file. Previous analysis had shown that, by increasing the release tolerance 
(reltol> 0.2), termination of the analysis could be avoided. However, the results had not been in good 
agreement with the benchmark [3,4], Therefore, additional stabilization had to be introduced in 
order to obtain agreement with the benchmark case. 
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The results computed when contact stabilization (cs) was added are plotted in Figure 40. A 
small stabilization factor ( 0 ?= lxl O' 6 ) was used for all cases since it had yielded good results in 
previous analyses [3,4], Initially, the release tolerance value was set at the default value (reltol= 0.2) 
(solid green line). For this parameter combination, the load increased up to the critical point, and 
delamination propagation started while the load remained constant (solid green line). Also for the 
stable propagation path, the result was in good agreement with the benchmark result (solid grey 
line). Then, the release tolerance was reduced to reltol= 0.1 (solid blue line) and further to 
reltol= 0.01 (solid red line). For all cases, the results were in good agreement with the benchmark 
results over the entire load/displacement range. 

An alternative way to plot the benchmark is shown in Figure 41 where the applied load P is 
plotted versus the increase in delamination length a*. This way of presenting the results may be of 
advantage for large structures where local delamination propagation may have little effect on the 
global stiffness of the structure and may therefore not be visible in a global load/displacement plot. 
However, extracting the delamination length a from the finite element results required more 
manual, time consuming, post-processing of the results compared to the relatively simple and 
readily available output of nodal displacements and forces. The results plotted in Figure 41 are the 
same as those that were discussed above and were shown as global load/displacement plots in 
Figure 40. The conclusions that can be drawn from the plots in Figure 41 are identical to those 
discussed in the above for Figure 40. 

A result obtained from a three-dimensional model is shown in Figures 42. Based on the results 
from two-dimensional planar models shown above, contact stabilization was added to help 
overcome convergence issues. A release tolerance value (reltol=0A) was used which had yielded 
good results previously. The load/displacement result computed when a small stabilization factor 
(cs=lxl0' 6 ) was added is plotted in Figure 42. The load increase stopped and delamination 
propagation started (solid blue line) shortly before reaching the critical point of the benchmark 
solution (solid grey line). As mentioned earlier, deviation from the benchmark may be explained by 
the fact that the benchmark results were created using two-dimensional planar finite element models 
of the MMB specimen. As shown in Figure 29, the three-dimensional models yield a failure index 
which is somewhat higher across the entire width than the results from two-dimensional planar 
finite element models. Therefore, delamination propagation is expected to start prior to the 
benchmark results obtained from two-dimensional planar finite element analysis, which ultimately 
leads to a shift of the entire results plot towards lower loads. 

5.4 Results from automated delamination propagation analysis for 80% mode II 

5.4.1 Computed delamination propagation for applied displacement 

As discussed above, the propagation analysis was performed in two steps using the models 
shown in Figures 4b and 6b starting from an initial delamination length, a<j=25A mm. In the first 
step, a displacement of w=1.5 mm was applied which nearly equaled the critical displacement, 
w C rit= 1.65 mm, determined earlier. In the second step, the applied displacement was increased to 
w=8.0 mm. For this second step, automatic incrementation was used in ABAQUS® and a small 

increment size (0.5% of the total step) was chosen at the beginning of the step. To reduce the risk of 
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numerical instability and early termination of the analysis, a minimum allowed increment size (10' 
of the total step) was also chosen. The analysis was limited to 5000 increments. 

Initially, analyses were performed using two-dimensional planar models (shown in Figure 4a) 
without stabilization or viscous regularization. Using the default value reltol= 0.2, the load and 
displacement exceeded the critical point and the analysis terminated (solid green line) as shown in 
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Figure 43 where the computed resultant force (load P) is plotted versus the applied displacement w. 
By increasing the release tolerance to reltol= 0.5, - as suggested by the ABAQUS® error in the 
message (.msg) file - it was possible to extend the analysis without an error message until the load 
drop occurred and delamination propagation started as shown in Figure 43, (solid light blue line). 
However, the load drop occurred past the critical point and the analysis terminated soon afterwards 
due to convergence problems. The release tolerance was not increased any further as suggested in 
the ABAQUS® error in the message (.msg) file. Previous analysis had shown that by increasing the 
release tolerance (reltol> 0.2), termination of the analysis could be avoided, however, the results had 
not been in good agreement with the benchmark [3,4], Therefore, it was decided to introduce 
additional stabilization to obtain better agreement with the benchmark case. 

The results computed when contact stabilization (cs) was added are plotted in Figure 44. A 
small stabilization factor (cv=lxlO" 6 ) was used for all cases since it had yielded good results in 
previous analyses [3,4], Initially, the release tolerance value was set at the default value (reltol= 0.2) 
(solid green line). For this parameter combination, the load increased past the critical point and 
delamination propagation started while the load dropped (solid green line). The load/displacement 
path then ran parallel to the constant displacement branch of the benchmark result and later 
followed the stable propagation branch of the benchmark result. To reduce the observed overshoot, 
the release tolerance was reduced. For a stabilization factor of cs=lxl 0" 6 and release tolerance 
reltol= 0.1 (solid blue line), the overshoot was reduced, and the computed load/displacement path 
then ran closer to the constant displacement branch of the benchmark result. However, just before 
reaching the minimum load, the displacement increased slightly followed by a small load drop. The 
results then followed the stable propagation branch of the benchmark result. Further reducing the 
release tolerance (reltoh 0.01) yielded results that were in excellent agreement with the benchmark 
(solid red line). Due to the fine mesh, only a small saw-tooth pattern was observed along the stable 
path. 

An alternative way to plot the benchmark is shown in Figure 45 where the applied displacement 
w is plotted versus the increase in delamination length a*. The results plotted in Figure 45 are the 
examples that were discussed above and were shown in the global load/displacement plots of 
Figure 44. The conclusions that can be drawn from this plot are identical to those discussed above. 

The result obtained from a three-dimensional analysis is shown in Figure 46. Based on the 
results from two-dimensional planar models shown above, contact stabilization was added to the 
analysis to help overcome convergence issues. Since the computer time increases with tighter 
release tolerances, i.e. smaller release tolerance values, first a release tolerance value (reltol= 0.1) 
was selected. The result (solid blue line) computed for a release tolerance value {reltol= 0.1) and a 
small stabilization factor (cv=lxl0‘ A ) are plotted in Figure 46. The load dropped and delamination 
propagation started upon reaching the critical point of the benchmark solution (solid grey circles and 
solid grey line). The load/displacement path then ran parallel to the benchmark result but was 
shifted slightly towards lower loads. The analysis was stopped after 6.81 10 5 CPU seconds (7.9 
days). To improve the results, the release tolerance value was set to reltol= 0.01. The results (solid 
red line) are practically identical to the results obtained for reltol= 0.01, however the computation 
time increased and the analysis was stopped after 1.97 1 0 6 CPU seconds (22.8 days). Deviation 
from the benchmark may be explained by the fact that the benchmark results were created using 
two-dimensional planar finite element models of the MMB specimen. As shown in Figure 29, the 
three-dimensional models yield a failure index, which is somewhat higher across the entire width 
than the results from two-dimensional planar finite element models. Therefore, delamination 
propagation is expected to start prior to the benchmark results obtained from two-dimensional 
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planar finite element analysis, which ultimately leads to a shift of the entire results plot towards 
lower loads. The current results for 80% mode II resemble the results obtained from the analysis of 
the ENF specimen (pure mode II) [4], 

5.4.2 Computed delamination propagation for applied quasi-static load 

The propagation analysis was performed in two steps using the models shown in Figures 4b and 
6b starting from an initial delamination length, ao=25.4 mm. In the first step, a load P=73 0 N was 
applied which equaled nearly the critical load, P cn r= 751 N, determined earlier. In the second step, 
the total load was increased (P=\ 000 N). Automatic incrementation was used with a small 
increment size at the beginning (0.5% of the total step) and a very small minimum allowed 
increment (10' of the total step) to reduce the risk of numerical instability and early termination of 
the analysis. The analysis was limited to 5000 increments. 

The same steps discussed in the section on applied displacement were followed. Initially, 
analyses were performed using two-dimensional planar models without stabilization or viscous 
regularization. Using the default value (reltol= 0.2), the load increase stopped after reaching the 
critical point, but the analysis terminated immediately (solid green line) due to convergence 
problems, as shown in Figure 47. The release tolerance was not increased as suggested by the 
ABAQUS® error in the message (.msg) file. Previous analysis had shown that by increasing the 
release tolerance (reltol> 0.2), termination of the analysis could be avoided, however, the results had 
not been in good agreement with the benchmark [3,4], Therefore, it was decided to introduce 
additional stabilization to obtain better agreement with the benchmark case. 

The results computed when contact stabilization (cs) was added are plotted in Figure 48. A 
small stabilization factor ( 0 ?= lxl O' 6 ) was used for all cases since it had yielded good results in 
previous analyses [3, 4], Initially, the release tolerance value was set at the default value (reltol= 0.2) 
(solid green line). For this parameter combination, the load increased up to the critical point and 
delamination propagation started while the load remained constant (solid green line). Also for the 
stable propagation path, the result was in good agreement with the benchmark result (solid grey 
line). Then, the release tolerance was reduced to reltol= 0.1 (solid blue line) and further to 
reltol= 0.01 (solid red line). For all cases, the results were in good agreement with the benchmark 
results over the entire load/displacement range. 

An alternative way to plot the benchmark is shown in Figure 49 where the applied load P is 
plotted versus the increase in delamination length a*. This way of presenting the results may be of 
advantage for large structures where local delamination propagation may have little effect on the 
global stiffness of the structure and may therefore not be visible in a global load/displacement plot. 
The results plotted in Figure 49 are the same as those that were discussed above and were shown as 
global load/displacement plots in Figure 48. The conclusions that can be drawn from the plots in 
Figure 49 are identical to those discussed in the above for Figure 48. 

A result obtained from a three-dimensional model is shown in Figures 50. Based on the results 
from two-dimensional planar models shown above, contact stabilization was added to help 
overcome convergence issues. Since the computer time increases with tighter release tolerances, i.e. 
smaller release tolerance values, first a release tolerance value reltol= 0. 1 was selected. The result 
(solid blue line) computed for a release tolerance value (reltol= 0.1) and a small stabilization factor 
(cs=lxl0' 6 ) are plotted in Figure 50. The load increase stopped (solid blue line) and delamination 
propagation started shortly before reaching the critical point of the benchmark solution (solid grey 
line). The load/displacement path then ran parallel to the benchmark result but was shifted slightly 
towards lower loads. To improve the results, the release tolerance value was set to reltol= 0.01. The 
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results (solid red line) are practically identical to the results obtained for reltol= 0.1, however the 
computational time increased. As discussed earlier, deviation from the benchmark may be explained 
by the fact that the benchmark results were created using two-dimensional planar finite element 
models of the MMB specimen. As shown in Figure 29, the three-dimensional models yield a failure 
index which is somewhat higher across the entire width than the results from two-dimensional 
planar finite element models. Therefore, delamination propagation is expected to start prior to the 
benchmark results obtained from two-dimensional planar finite element analysis, which ultimately 
leads to a shift of the entire results plot towards lower loads. The current results for 80% mode II 
resemble the results obtained from the analysis of the ENF specimen (pure mode II) [4], 

The results from above were plotted for three distinct locations across the width of the specimen 
(two edges and center as shown later in Figure 56) as shown in Figure 51, where the applied load P 
is plotted versus the increase in delamination length a*. This way of presenting the results is shown 
since it may be of advantage for large structures where local delamination propagation may have 
little effect on the global stiffness of the structure, and may therefore not be visible in a global 
load/displacement plot. However, extracting the delamination length a from finite element results 
required more manual, time consuming post-processing of the results - especially for the three- 
dimensional models - compared to the relatively simple and readily available output of nodal 
displacements and forces. 

Initially during the unstable growth phase, the results for both simulations (reltol=0. 1 , cs=lxl 0‘ 6 
solid lines and reltol= 0.01, cs= 1 x 1 O' 6 dashed lines) and all three locations (green line for edge 1, 
blue line for center location and red line for edge 2) were practically identical. All the result curves 
ran parallel to the benchmark result (solid grey line) but were shifted slightly towards lower loads. 
The unstable phase, however created longer delaminations compared to the benchmark result. For 
the stable path of delamination propagation, small differences can be observed between the two 
simulations (solid versus dashed lines). In this respect, more detailed information about the 
propagation can be obtained from the plot where the applied load is plotted versus the increase in 
delamination length than from the traditional load/displacement plot. The current results for 80% 
mode II resemble the results obtained from the analysis of the ENF specimen (pure mode II) [4], 

5.5 Computed delamination front shape 

Besides matching the load displacement behavior of benchmark results, a delamination 
propagation analysis should also yield a delamination front shape that is representative of the actual 
failure. An example of delamination front shapes observed by opening tested MMB specimens are 
shown in Figure 52 [7]. From the initial straight delamination front, which is formed by the edge 
of the polytefrafluoroethylene (PTFE) insert, the delamination develops into a curved thumbnail 
shaped front (blue line) as shown in Figure 52a for 20% mode II. The front is similar to fronts 
observed in the mode I DCB specimen [3]. A pronounced thumbnail shaped front is observed for 
the pre-cracks (red line) shown in Figures 52 b and c for 50% and 80% mode II, respectively. The 
final front (blue line) exhibits less curvature as would be expected from the increased mode II 
contribution. The fronts are somewhat jagged, suggesting that growth happens in one location then 
stops and continues at another location across the width. On average, however, the delaminations 
appear to grow uniformly across the width. 

A straight front across almost the entire width of the specimen is to be expected when looking at 
the failure index distribution plotted in Figure 29. However, due to the drop of the failure index near 
the edges, delamination propagation is expected to locally lag behind and turn the initially straight 
front into a thumbnail front for Gh/Gt=0.2 and 0.5 similar to fronts observed in a Double Cantilever 
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Beam (DCB) specimen. For Gu/Gt =0.8, the drop of the failure index near the edges is less 
pronounced, and a straight front across the entire width is expected. 

Delamination propagation computed using the model with a uniform mesh across the width 
(Figure 5) is shown in Figure 53 for Gu/Gt =0.2. Plotted on the bottom surface (defined in Figure 5) 
are the contours of the bond state, where the delaminated section appears in blue and the intact 
(bonded) section in red. The transition between the colors (in green/yellow) indicates the location of 
the delamination front. The initial straight front is shown in Figure 53a and was added for 
clarification in Figures 53d and e. The first propagation was observed towards the edges of the 
specimen (Figure 53b) as expected from the distribution of the failure index (Figure 29) where two 
local maxima are observed at y/B ± 0.33. Then, as expected, the front advanced more in the center 
and lagged behind at the edges as shown in Figure 53c. Later in the analysis, the delamination 
continued to propagate through the entire specimen basically as a straight front as shown in 
Figures 53d and 53e. Locally however, always one element at each edge lagged behind as shown in 
Figures 53c and 53d. After the shape stabilized, new propagation usually started in the center first as 
shown in Figure 53e, the exact location being determined by the maximum of the failure index 
distribution across the width. 

Delamination fronts computed using the model with a uniform mesh across the width (Figure 5) 
are shown in Figures 54 and 55 for Gu/Gt =0.5. For the case of applied displacement, w, the 
propagation was observed first in the center section of the specimen (Figure 54a) rather than at the 
location of the two local maxima in the failure index distribution (Figure 29) observed at y/B ± 0.33. 
Then, the front advanced across the entire width but lagged behind at the edges as shown in 
Figure 54b. Later, the delamination propagated across the width of the specimen as a straight front 
shape with only one element at each edge lagging behind as shown in Figure 54c. The final front as 
shown in Figure 54d, appeared jagged locally but on average moved forward as a straight front 
across the entire width, which may be explained by the increase in mode II for longer delamination 
lengths (increased Gu/Gt as shown in Figure 25). Delaminations under pure mode II generally 
exhibit straight fronts [4,10]. For the case of applied load, P, the propagation was observed first in 
the center section of the specimen (Figure 55a). Then, the front advanced across the entire width but 
lagged behind at the edges as shown in Figure 55b. Later, the delamination appeared jagged locally 
but on average moved forward as a straight front across the entire width as shown in Figures 55c 
and 55d, which may again be explained by the increase in mode II for longer delamination lengths 
(increased Gu/Gt as shown in Figure 25). 

Delamination propagation computed using the model with a uniform mesh across the width 
(Figure 5) is shown in Figure 56 for Gu/Gt =0.8. The initial straight front is shown in Figure 56a 
and 56b and was added for clarification in Figures 56d and 56e. The first propagation was observed 
towards the edges of the specimen (Figure 56c) as expected from the distribution of the failure 
index (Figure 29) where two local maxima are observed at y/B ±0.33. Then, the delamination 
appeared jagged locally but on average moved forward as a straight front across the entire width as 
shown in Figures 56d, which may be explained by the already high mode II ratio ( Gi/Gt =0.8) for 
this configuration. The front remained jagged until the end of the analysis as shown in Figure 56e. 
The current results for 80% mode II resemble the results obtained from the analysis of the ENF 
specimen (pure mode II) [4], 

Computed delamination fronts obtained from three-dimensional solid models generally matched 
the experimentally observed shapes shown in Figure 52. Improved results may be obtained by 
refining the mesh across the width. 


16 


6. SUMMARY AND CONCLUSIONS 

The development of benchmark examples, which allow the assessment of the static 
delamination propagation prediction capabilities was presented and demonstrated for ABAQUS® 
Standard. The examples are based on finite element models of the mixed-mode I/II Mixed-Mode 
Bending (MMB) specimen. The models are independent of the analysis software used and allow the 
assessment of the automated delamination growth prediction capabilities in commercial finite 
element codes based on the virtual crack closure technique (YCCT). 

First, quasi-static benchmark results were created based on the approach developed in reference 
3. Two-dimensional finite element models were used for simulating the MMB specimens with 
different mixed-mode ratio (20%, 50% and 80% mode II) and different initial delamination lengths, 
ao. For each delamination length modeled, the load and displacements were monitored. The mixed- 
mode I/II strain energy release rate was calculated for a fixed applied displacement. It was assumed 
that the delamination propagated when the total strain energy release rate reached the fracture 
toughness value. Thus, critical loads and critical displacements for delamination propagation were 
calculated for each initial delamination length modeled. From these critical load/displacement 
results, benchmark solutions were created. It was assumed that the load/displacement relationship 
computed during automatic propagation should closely match the benchmark cases. 

After creating the benchmark cases, the approach was demonstrated for the commercial finite 
element code ABAQUS®. Starting from an initially straight front, the delamination was allowed to 
propagate under quasi-static loading based on the algorithms implemented into the software. Input 
control parameters were varied to study the effect on the computed delamination propagation and 
growth. 

The results showed the following: 

• The benchmarking procedure proved valuable by highlighting the issues associated 
with choosing the input parameters of the particular implementation. 

• In general, good agreement between the results obtained from the automated 
propagation analysis and the benchmark results could be achieved by selecting input 
parameters that had previously been determined during analyses of mode I Double 
Cantilever Beam and mode II End Notched Flexure specimens. 

• In particular, the results for automated delamination propagation analysis under quasi- 
static loading showed the following: 

o Using the default release tolerance (reltol= 0.2) as suggested in the ABAQUS® 
handbook or increasing the value, as suggested in the user's manual, may help 
to overcome convergence problems, however, it leads to an undesired 
overshoot of the computed result compared to the benchmark. 

o A combination of release tolerance and contact stabilization is required to 
obtain more accurate results. 

o A gradual reduction of the release tolerance and contact stabilization over 
several analyses is suggested. 

o Good agreement between analysis results and the benchmarks could be achieved 
for release tolerance values (reltol< 0.1) in combination with contact stabilization 
(cs=lxl0" 6 ). 

o Computed delamination fronts obtained from three-dimensional solid models 
generally matched the experimentally observed shapes. Improved results may be 
obtained by refining the mesh across the width. 
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o The analyses for three-dimensional models of a simple MMB specimen required 
several days of computation time. Improvements to the implementation are 
required to make analysis of larger structural components computationally 
affordable. 

Overall, the results are promising and the current findings concur with previously published 
conclusions [3,4,5]. However, further assessment for mixed-mode delamination fatigue onset and 
growth is required. Additional studies should also include the assessment of the propagation 
capabilities in more complex mixed-mode specimens and on a structural level. 

Assessing the implementation in one particular finite element code illustrated the value of 
establishing benchmark solutions since each code requires specific input parameters unique to its 
implementation. Once the parameters have been identified, they may then be used as a starting point 
to model delamination growth for more complex configurations. 
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TABLE I. MATERIAL PROPERTIES [4]. 


Unidirectional Graphite/Epoxy Prepreg 

£ll = 161 GPa 

£22 = 11.38 GPa 

£33 = 11.38 GPa 

V12 = 0.32 

V13 = 0.32 

V23 = 0.45 

G12 = 5.2 GPa 

G13 = 5.2 GPa 

G23 = 3.9 GPa 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes the ply principal 
axis that coincides with the direction of maximum in-plane Young’s modulus (fiber direction). Index 22 denotes the 
direction transverse to the fiber in the plane of the lamina and index 33 the direction perpendicular to the plane of the 
lamina. 


TABLE II. FRACTURE PARAMETERS. 


Fracture Toughness Data [7,10,1 1] - Figures 2, A1 

Gj c = 0.212 kJ/m 2 Gjj c = Gjji c =0.774 kJ/m 2 rj= 2. 1 
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APPENDIX 

Delamination propagation analysis in ABAQUS® 

Delamination propagation at the interfaces in laminated composites can be simulated in 
ABAQUS®[14], The interface along which the delamination (or crack) propagates must be 
indicated in the model using a fracture criterion definition. The fracture energy release rates at 
the crack tips in the interface elements are calculated based on the virtual crack closure technique 
(VCCT) [14,15], 

Required input for ABAQUS® 

The input required to perform a delamination propagation analysis under quasi-static load in 
ABAQUS® Standard is discussed in the following paragraphs. It is assumed that the reader is 
familiar with ABAQUS® Standard and the syntax used in the input file (.inp). The focus is 
therefore on the specific input that relates to delamination propagation in ABAQUS®. An 
example input file is given at the end of this appendix to provide an overview of an entire 
analysis and assist the readers in creating their own analyses. 

Input for delamination propagation 

The interface along which the delamination (or crack) propagates must be indicated in the 
model using a fracture criterion definition: 

* DEBOND , SLAVE=VCCT_TOP , MASTER=VCCT_BOT , FREQ=1 

♦FRACTURE CRITERION, TYPE=VCCT , MIXED MODE BEHAVIOR=BK, TOLERANCE=<reltol> 

<GI c> , <GI I c> , <GI 1 1 c> , <e ta> 

* CONTACT CONTROLS , STABILIZE, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
l.E-06, 0, 0.1 

where vcct_top and vcct_bot are interface surfaces as shown in Figure Al, and <reitoi> is 
the release tolerance within which the crack propagation criterion must be satisfied. The critical 
energy release rates <gic>,<giic>,<giiic> and the curve fit parameter <eta> are obtained from 
the mixed-mode failure criterion as shown in Figure A2. A quasi-static mixed-mode fracture 
criterion is determined for a material by plotting the interlaminar fracture toughness, G c , versus the 
mixed-mode ratio, Gji/Gt as discussed earlier and shown in Figure A2. The input for contact 
stabilization is provided in the contact controls definition. 

Example input files 

An input file is given to provide an overview of an entire analysis and assist the readers in 
creating their own analyses. The analysis was divided into two steps. In the first step, a 
displacement just below the critical value, w crih was applied to avoid delamination propagation in 
this step. Thus, large increments could be used to get up to the critical point. In the second step, 
the final desired displacement was defined. Automatic incrementation was used with a small 
increment size at the beginning (0.5% of the total step) and a very small minimum allowed 
increment (10' of the total step) to reduce the risk of numerical instability and early termination of 
the analysis. The analysis was limited to 5000 increments. 

For all analyses, the input to define the fracture criterion (<Gic>,<Giic>,<Giiic>,<eta>), 
were kept constant. The ABAQUS® keywords shown in bold type were discussed in detail in the 
previous paragraphs. 
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Input file for fatigue onset and growth analysis 

*HEADING 

MMB-UD-IM7/8552, a=30 mm, ratio=0.2 
units: mm, N, MPa 

*** elements , nodes, material , etc 
*ELEMENT, TYPE=CPE4I , ELSET=EALL 


*NSET, NSET=BONDED, GENERATE 
1069, 1533, 8 

* * 

**** surface and contact definition **** 

* SURFACE, TYPE=ELEMENT , NAME=VCCT_BOT 
EL_BOT , S3 

* SURFACE, TYPE=ELEMENT , NAME=VCCT_TOP 
EL_TOP, SI 

*INITIAL CONDITION, TYPE=CONTACT 
VCCT_TOP, VCCT_BOT, BONDED 

*CONTACT PAIR, INTERACTION=VCCT, ADJUST=BONDED 
VCCT_TOP, VCCT_BOT 

-k * * * VCC T STAGE I **** 

* SURFACE INTERACTION, NAME=VCCT 
<width> 

**** rigid elements to simulate MMB fixture 
*** extra nodes for rigid elements 
*NSET, NSET=MMBREF, UNSORTED 
**** mmb fixture loading 

100006, 

* * 

^NSET , NSET=NPIN, UNSORTED 
**** MMB rollers 

100001 , 

* * 

^NODE , NSET=NFIXTURE 

100001, 2.25, 152.399994, 0. 

100002, 34., 152.399994, 0. 

100003, 2.25, 102.000168, 0. 

100004, 34., 102.000168, 0. 

** ratio=0.2 

100005, 34., 9.05, 0., 

100006, 3.25, 9.05, 0., 

** ratio=0.5 

**100005, 34., 60.72, 0., 

**100006, 3.25, 60.72, 0., 

** ratio=0 . 8 

**100005, 34., 74.75, 0., 

**100006, 3.25, 74.75, 0., 

*** rigid beam elements 

*ELEMENT , TYPE=R2D2 , ELSET=MMB_RIG 

100001 , 8 , 100002 

100002, 100002, 100004 

**100003, 272, 100004 

100003, 100003, 100004 

100004, 100004, 100005 

100005, 100005, 100006 
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FRIGID BODY, ELSET=MMB_RIG, REF N 0 DE =MMBRE F , PIN NSET=NPIN 
* * 

** mode II loading 
^EQUATION 
2 

616,1,1.0,100003,1,-1.0 

**** VCCT input 
^parameter 

** fracture toughness: 

GIc = 0.212 
GIIc = 0.774 
GIIIc = 0.774 
** B-K parameter: 
eta=2 . 1 

** Damage and tolerance parameters 
damv=0 . 0 
reltol=0 . 01 

** width in the plane stress/strain direction 
width =25.4 

* * * 

********************** STEP 1 ******************************** 

*STEP, NLGEOM, INC=500 

* STATIC 

0.1, 1., 0.000000001, 1. 

**** MMB loading 

*** applied displacement w for MMB, 2D- full model 
* BOUNDARY 

100006, 1, 1, -1.5 

* * 

*CONTROLS , PARAMETERS=TIME INCREMENTATION 

r r i r r r r ^ 0 

*** fracture analysis using VCCT 

* DEBOND , SLAVE=VCCT_TOP ,MASTER=VCCT_BOT , FREQ=1 

* FRACTURE CRITERION , TYPE=VCCT , MIXED MODE BEHAVIOR=BK , TOLERANCE=<reltol> 
<GI c> , <GI I c> , <GI 1 1 c> , <e ta> 

** 

***CONTACT CONTROLS, STABILIZE, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
** l.E-06, 0, 0.1 

* * 

^OUTPUT, FIELD, VARIABLE=PRESELECT , FREQ=1 
^CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
dbt, dbsf , dbs, openbc, crsts, enrrt, ef enrrtr , bdstat 
* * 

^Output, history, VARIABLE=PRESELECT, freq=l 
*NODE output, NSET=MMBREF 
RF1 

*NODE output, NSET=MMBREF 
U1 

***NODE output, NSET=NSUP 
** RF1 , RF2 

***NODE output, NSET=NSUP 
** Ul, U2 

^CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP, NSET=BONDED 
enrrt 

*END STEP 

********************** STEP 2 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

*STEP, NLGEOM, INC= 5000 
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* STATIC 

0.00500000 , 1.0000 , 1 . 0000000E-18 , 1.0 

**** MMB loading 

*** applied displacement w for MMB, 2D- full model 
* BOUNDARY 

100006 , 1, 1, -8.0 

* * 

♦CONTROLS , PARAMETERS=TIME INCREMENTATION 

r r r r r r r ^ 0 

*CONTACT PRINT 

*** fracture analysis using VCCT 

* * *DEBOND , SLAVE=VCCT_TOP , MASTER=VCCT_BOT , FREQ=1 
*DEBOND , SLAVE=VCCT_TOP ,MASTER=VCCT_BOT , FREQ=1 

* FRACTURE CRITERI ON, TYPE=VCCT, MIXED MODE BEHAVIOR=BK, TOLERANCE=<reltol> 
<GIc> , <GI Ic> , <GII Ic> , <eta> 

♦CONTACT CONTROLS, STABILIZE, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
l.E-06, 0, 0.1 

* * 

♦OUTPUT, FIELD, VARIABLE=PRESELECT, FREQ=1 
♦CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP 
dbt, dbsf , dbs, openbc, crsts, enrrt, efenrrtr, bdstat 
♦Output, history, VARIABLE=PRESELECT, freq=l 
♦NODE output, NSET=MMBREF 
RF1 

♦NODE output, NSET=MMBREF 
U1 

♦CONTACT OUTPUT, MASTER=VCCT_BOT , SLAVE=VCCT_TOP, NSET=BONDED 
enrrt 

*END STEP 
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dimensions 

b 25.4 mm 

2h 4.5 mm 

2L 100.8 mm 

2L* 152.4 mm 

a Q 25.4 mm 

layup: IM7/8552 [0] 24 


Figure 1. Mixed Mode Bending specimen (MMB) ( dimensions from [7]). 
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Figure 2. Mixed mode fracture criterion for IM7/8552. 
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(a). View of full model - G ^/Gj-0.2, c=92.9 mm. 

rigid elements to 



z,w 

A 


multi-point constraints 
to enforce same w 
displacement and allow u 
displacement (sliding) 


x,u 



rigid elements to 
simulate load apparatus 

bonded nodes 


(c). Detail of load introduction zone. 

Figure 3. Deformed 2D FE-model of a MMB specimen with mixed mode ratio G n /G T =0.2. 
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rigid elements to 
simulate load 
apparatus 



(a). Mixed mode ratio Gjj/G t = 0.5, c=41.3 mm. 


c 



(b). Mixed mode ratio GjfG T =0.8, c=27.3 mm. 


Figure 4. Deformed 2D FE models ofMMB specimens with different mixed-mode ratios. 
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delamination front 


bonded 

nodes 


top surface 
bottom surface 


(a). View of full model - G n /G T =0.2, c=92.9 mm. 



(b). Detail of specimen tip and crack tip zone. 

Figure 5. Deformed 3D FE-model of a MMB specimen with mixed mode ratio G U /G T =0.2. 
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Figure 6. Deformed 3D FE-models of a MMB specimen with different mixed mode ratios. 
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applied displacement w, mm 

Figure 7. Load-displacement behavior for MMB specimens (20% mode 11) 
with different delamination lengths a^ 
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Figure 8. Use of mixed-mode fracture criterion to obtain G for different mode ratios G /G . 

c II T 
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Figure 9. Calculated critical load-displacement behavior for a MMB specimen (20% mode II). 




applied displacement w, mm 

Figure 10. Calculated critical behavior and resulting benchmark case for 
applied displacement, w, for a MMB specimen (20% mode II). 
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Figure 1 1 . Benchmark case for applied displacement, w, plotted versus 
increase in delamination length, a*, for a MMB specimen (20% mode II). 
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Figure 12. Load-displacement behavior for MMB specimens (50% mode II) 


with different delamination lengths a f . 
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Figure 13. Calculated critical load-displacement behavior for a MMB specimen (50% mode II). 
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Figure 14. Calculated critical behavior and resulting benchmark cases for 
applied load, P, and displacement, w,for a MMB specimen (50% mode II). 
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Figure 15. Benchmark case for applied displacement, w, plotted versus 
increase in delamination length, a*, for a MMB specimen (50% mode II). 


applied 
load 
P, N 



Figure 16. Benchmark case for applied load, P, plotted versus increase 
in delamination length, a*, for a MMB specimen (50% mode II). 
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Figure 17. Load-displacement behavior for MMB specimens (80% mode II) 
with different delamination lengths, a^. 


10.0 


load P, 



applied displacement w, mm 

Figure 18. Calculated critical load-displacement behavior for a MMB specimen (80% mode II). 
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Figure 19. Calculated critical behavior and resulting benchmark cases for 
applied load, P, and displacement, w, for a MMB specimen ( 80% mode II). 
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Figure 20. Benchmark case for applied displacement, w, plotted versus increase in 
delamination length, a*, for a MMB specimen (80% mode II). 


35 



Figure 21 . Benchmark case for applied load, P, plotted versus increase in 
delamination length, a*, for a MMB specimen (80% mode II). 
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Figure 22. Computed total energy release rate. 
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Figure 23. Computed mode II energy release rate component. 
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Figure 24. Computed mixed-mode ratio. 
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Figure 25. Dependence of computed mixed-mode ratio on delamination length. 
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Figure 26. Computed total energy release rate distribution across the width of a MMB specimen. 
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Figure 27. Computed mode II energy release rate distribution 
across the width of a MMB specimen. 
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Figure 28. Computed mixed-mode ratio distribution across the width of a MMB specimen. 
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Figure 29. Failure index distribution across the width of a MMB specimen. 
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Figure 30. Computed critical load-displacement behavior for a MMB specimen (20% mode II) 
obtained from two-dimensional planar models ( CPE4I) with different release tolerance settings. 
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applied displacement w, mm 

Figure 31. Computed critical load-displacement behavior for a MMB specimen (20% mode II) 
obtained from two-dimensional planar models ( CPE4I) with added contact stabilization. 
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Figure 32. Computed displacement-propagation behavior for a MMB specimen (20% mode II) 
obtained from two-dimensional planar models ( CPE4I) with added contact stabilization. 
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applied displacement w, mm 

Figure 33. Computed critical load-displacement behavior forMMB specimen (20% mode II) 
obtained from three-dimensional solid models ( C3D8I) with added contact stabilization. 



increase in delamination length a* mm 

Figure 34. Computed displacement-propagation behavior for a MMB specimen (20% mode II) 
obtained from three-dimensional solid models ( C3D8I) with added contact stabilization. 
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applied displacement w, mm 

Figure 35. Computed critical load-displacement behavior for MMB specimen (50% mode II) 
obtained from two-dimensional planar models ( CPE4I) with different release tolerance settings. 
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Figure 36. Computed critical load-displacement behavior for MMB specimen ( 50% mode II) 
obtained from two-dimensional planar models ( CPE4I) with added contact stabilization. 


43 


applied 
displacement 
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Figure 37. Computed critical displacement-propagation behavior for MMB specimen (50% mode 11) 

obtained from two-dimensional planar models. 



applied displacement w, mm 

Figure 38. Computed critical load-displacement behavior for MMB specimen (50% mode II) 
obtained from three-dimensional solid models ( C3D8I) with added contact stabilization. 
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Figure 39. Computed critical load-displacement behavior for MMB specimen (50% mode II) 
obtained from two-dimensional planar models ( CPE4I) with different release tolerance settings. 
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Figure 40. Computed critical load-displacement behavior for MMB specimen (50% mode II) 
obtained from two-dimensional planar models ( CPE4I) with added contact stabilization. 
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Figure 41. Computed critical load-propagation behavior for MMB specimen (50% mode II) 

obtained from two-dimensional planar models. 



displacement w, mm 

Figure 42. Computed critical load-displacement behavior for MMB specimen ( 50% mode II) 
obtained from three-dimensional solid models ( C3D8I) with added contact stabilization. 
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Figure 43. Computed critical load-displacement behavior for MMB specimen ( 80% mode II) 
obtained from two-dimensional planar models ( CPE4I) with different release tolerance settings. 
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Figure 44. Computed critical load-displacement behavior for MMB specimen ( 80% mode II) 
obtained from two-dimensional planar models ( CPE4I) with added contact stabilization. 
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Figure 45. Computed critical displacement-propagation behavior for MMB specimen (80% mode II) 

obtained from two-dimensional planar models. 
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Figure 46. Computed critical load-displacement behavior for MMB specimen ( 80% mode II) 
obtained from three-dimensional solid models ( C3D8I) with added contact stabilization. 
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Figure 47. Computed critical load-displacement behavior for MMB specimen ( 80% mode II) 
obtained from two-dimensional planar models ( CPE4I) with different release tolerance settings. 
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Figure 48. Computed critical load-displacement behavior for MMB specimen ( 80% mode II) 
obtained from two-dimensional planar models ( CPE4I) with added contact stabilization. 
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Figure 49. Computed critical load-propagation behavior for MMB specimen (80% mode II) 

obtained from two-dimensional planar models. 



displacement w, mm 

Figure 50. Computed critical load-displacement behavior for MMB specimen ( 80% mode II) 
obtained from three-dimensional solid models ( C3D8I) with added contact stabilization. 
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Figure 51. Computed critical load-propagation behavior for MMB specimen (80% mode II) 
obtained from three-dimensional solid models ( C3D8I) with added contact stabilization. 
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a. 20 % mode II. 



b. 50% mode II. 



c. 80% mode II. 


Figure 52. Fracture surfaces of typical MMB specimens [7J. 
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a. Initial delamination front shape (Bottom surface ofFE model in Figure 5). 



b. Detail of first delamination propagation. c. Detail of delamination propagation. 
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Figure 53. Computed delamination front shape for MMB specimen (20% mode II). 
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c. Detail of advanced propagation. d. Detail of final front. 

Figure 54. Computed delamination front for MMB specimen (50% G^- applied displacement ). 
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c. Detail of advanced propagation. d. Detail of final front. 

Figure 55. Computed delamination front for MMB specimen (50% G n - applied load). 
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a. Initial delamination front shape (Bottom surface ofFE model in Figure 5). 




b. Detail of first initial delamination front. 


c. Detail of first delamination propagation. 



d. Detail of delamination propagation. e. Detail of delamination propagation. 

Figure 56. Computed delamination front shape for MMB specimen (80% mode II). 
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(b). Detail of specimen tip and crack tip zone. 



(c). Detail of load introduction zone. 


bonded nodes 
(BONDED) 


Figure Al. Deformed 2D FE-model of a MMB specimen. 
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Figure A2. Mixed mode fracture criterion for IM7 18552. 
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